﻿using System;

namespace MathService.DSP.Fourier
{
    /// <summary>Быстрое преобразование Фурье</summary>
    public static class FFT
    {
        /// <summary>Прямое преобразование отсчётов функции в спектр</summary>
        /// <param name="Values">Массив отсчётов функции</param>
        public static Complex[] FastFurierTransform(this double[] Values)
        {
            var lv_N = Values.Length;
            if(!lv_N.IsPowerOf2())
                lv_N = 1 << ((int)Math.Log(lv_N, 2) + 1);

            var lv_DoubleSpectrum = new double[lv_N * 2];
            var i_VLength = Values.Length;
            for(var i = 0; i < i_VLength; i++)
                lv_DoubleSpectrum[2 * i] = Values[i];

            fft(ref lv_DoubleSpectrum, false);

            var lv_ComplexSpectrum = new Complex[Values.Length];
            var i_CLength = lv_ComplexSpectrum.Length;
            for(var i = 0; i < i_CLength; i++)
                lv_ComplexSpectrum[i] = new Complex(lv_DoubleSpectrum[2 * i], lv_DoubleSpectrum[2 * i + 1]);

            return lv_ComplexSpectrum;
        }

        /// <summary>Прямое преобразование отсчётов функции в спектр</summary>
        /// <param name="Values">Массив отсчётов функции</param>
        public static Complex[] FastFurierTransform(this Complex[] Values)
        {
            var lv_N = Values.Length;
            if(!lv_N.IsPowerOf2())
            {
                var lv_Log2 = Math.Log(lv_N, 2);
                lv_Log2 -= Math.Round(lv_Log2);
                lv_N += (int)Math.Pow(2, lv_Log2);
            }

            var lv_Spectrum = new double[lv_N * 2];
            var lv_I = Values.Length;
            for(var i = 0; i < lv_I; i++)
            {
                lv_Spectrum[2 * i] = Values[i].Re;
                lv_Spectrum[2 * i + 1] = Values[i].Im;
            }

            fft(ref lv_Spectrum, false);

            var Spectrum = new Complex[lv_N];
            lv_I = Spectrum.Length;
            for(var i = 0; i < lv_I; i++)
                Spectrum[i] = new Complex(lv_Spectrum[2 * i], lv_Spectrum[2 * i + 1]);


            return Spectrum;
        }

        /// <summary>Обратное преобразование отсчётов спектра в отсчёты сигнала</summary>
        /// <param name="Spectrum">Массив отсчётов спектра</param>
        public static Complex[] FastFurierInverse(this Complex[] Spectrum)
        {
            var lv_N = Spectrum.Length;
            if(!lv_N.IsPowerOf2())
            {
                var lv_Log2 = Math.Log(lv_N, 2);
                lv_Log2 -= Math.Round(lv_Log2);
                lv_N += (int)Math.Pow(2, lv_Log2);
            }

            var lv_DoubleValues = new double[lv_N * 2];
            var lv_I = Spectrum.Length;
            for(var i = 0; i < lv_I; i++)
            {
                lv_DoubleValues[2 * i] = Spectrum[i].Re;
                lv_DoubleValues[2 * i + 1] = Spectrum[i].Im;
            }

            fft(ref lv_DoubleValues, true);

            var lv_ComplexValues = new Complex[lv_N];
            lv_I = lv_ComplexValues.Length;
            for(var i = 0; i < lv_I; i++)
                lv_ComplexValues[i] = new Complex(lv_DoubleValues[2 * i], lv_DoubleValues[2 * i + 1]);


            return lv_ComplexValues;
        }

        private static void fft(ref double[] Values, bool IsInverse)
        {
            var lv_N = Values.Length / 2;
            if(!lv_N.IsPowerOf2())
                throw new ArgumentException("Число элементов выборки должно быть степенью двойки");

            // ReSharper disable TooWideLocalVariableScope
            // ReSharper disable JoinDeclarationAndInitializer
            int jj;
            int m_max;
            int m;
            int i_step;
            int i;
            double w_temp;
            double w_r;
            double w_pr;
            double w_pi;
            double w_i;
            double theta;
            double temp_r;
            double temp_i;
            // ReSharper restore JoinDeclarationAndInitializer
            // ReSharper restore TooWideLocalVariableScope

            var i_sign = IsInverse ? -1 : 1;

            var n = lv_N << 1;
            var j = 1;


            int ii;
            //Операция бабочка
            for(ii = 1; ii <= lv_N; ii++)
            {
                i = (ii << 1) - 1;
                if(j > i)
                {
                    temp_r = Values[j - 1];
                    temp_i = Values[j];
                    Values[j - 1] = Values[i - 1];
                    Values[j] = Values[i];
                    Values[i - 1] = temp_r;
                    Values[i] = temp_i;
                }
                m = n >> 1;
                while(m >= 2 && j > m)
                {
                    j -= m;
                    m >>= 1;
                }
                j += m;
            }

            m_max = 2;
            var thetta0 = Consts.pi2 * i_sign;
            while(n > m_max)
            {
                i_step = m_max << 1;
                theta = thetta0 / m_max;
                //theta = Consts.pi2 / (i_sign * m_max);
                w_pr = Math.Sin(.5 * theta);
                w_pr *= -2 * w_pr;
                w_pi = Math.Sin(theta);
                w_r = 1;
                w_i = 0;
                for(ii = 1; ii <= (m_max >> 1); ii++)
                {
                    m = (ii << 1) - 1;
                    for(jj = 0; jj <= (n - m) / i_step; jj++)
                    {
                        i = m + jj * i_step;
                        j = i + m_max;
                        temp_r = w_r * Values[j - 1] - w_i * Values[j];
                        temp_i = w_r * Values[j] + w_i * Values[j - 1];
                        Values[j - 1] = Values[i - 1] - temp_r;
                        Values[j] = Values[i] - temp_i;
                        Values[i - 1] = Values[i - 1] + temp_r;
                        Values[i] = Values[i] + temp_i;
                    }
                    w_temp = w_r;
                    w_r = w_r * w_pr - w_i * w_pi + w_r;
                    w_i = w_i * w_pr + w_temp * w_pi + w_i;
                }
                m_max = i_step;
            }

            if(!IsInverse)
                for(i = 1; i <= 2 * lv_N; i++)
                    Values[i - 1] = Values[i - 1] / lv_N;
        }

        /// <summary>
        /// Целочисленное преобразование Фурье
        /// </summary>
        /// <param name="a">Массив целых чисел</param>
        /// <param name="invert">Обратное преобразование</param>
        public static void FFT_int(this int[] a, bool invert = false)
        {
            var n = a.Length;

            for(int i = 1, j = 0; i < n; i++)
            {
                var bit = n >> 1;
                for(; j >= bit; bit >>= 1)
                    j -= bit;

                j += bit;

                if(i >= j) continue;
                a[i] ^= a[j];
                a[j] ^= a[i];
                a[i] ^= a[j];
            }

            const int root = 5;
            const int root_1 = 4404020;
            for(var len = 2; len <= n; len <<= 1)
            {
                var wlen = invert ? root_1 : root;

                const int root_pw = 1 << 20;
                const int mod = 7340033;

                for(var i = len; i < root_pw; i <<= 1)
                    wlen = (int)(wlen * 1 * wlen % mod);
                for(var i = 0; i < n; i += len)
                {
                    var w = 1;
                    for(var j = 0; j < len / 2; ++j)
                    {
                        var u = a[i + j];
                        var v = (int)(a[i + j + len >> 1] * 1 * w % mod);
                        a[i + j] = u + v < mod ? u + v : u + v - mod;
                        a[i + j + len / 2] = u - v >= 0 ? u - v : u - v + mod;
                        w = (int)(w * 1 * wlen % mod);

                    }
                }

                if(invert)
                {
                    var nrev = reverse(n, mod);
                    for(var i = 0; i < n; i++)
                        a[i] = (int)(a[i] * 1 * nrev % mod);
                }

            }
        }

        private static int reverse(int n, int Mod)
        {
            throw new NotImplementedException();
        }
    }
}
